Quantum turbulence in condensate collisions: an application of the classical field method 
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We apply the classical field method to simulate the production of correlated atoms during the collision of two 
Bose-Einstein condensates. Our non-perturbative method includes the effect of quantum noise, and provides 
for the first time a theoretical description of collisions of high density condensates with very large out-scattered 
fractions. Quantum correlation functions for the scattered atoms are calculated from a single simulation, and 
show that the correlation between pairs of atoms of opposite momentum is rather small. We also predict the ex- 
istence of quantum turbulence in the field of the scattered atoms — a property which should be straightforwardly 
measurable. 
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In the same way as a classical electromagnetic field obey- 
ing Maxwell's equations arises as an assembly of photons all 
in the same quantum state, a Bose-Einstein condensate, com- 
posed of Bosonic atoms all in the same quantum state, be- 
haves very much like a classical field ^(x, t), whose equation 
of motion is the Gross-Pitaevskii equation 
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Nevertheless, there are phenomena in which the quantized na- 
ture of this field is important — for example, when two Bose- 
Einstein condensates collide at a sufficiently high velocity, a 
halo of elastically scattered atoms is produced 1 1, 2, 3J. The 
Gross-Pitaevskii equation with initial conditions correspond- 
ing to two Bose-Einstein condensates does not predict this 
scattering — it is a direct effect of the fact that the quantized 
field consists of interacting particles. 

Theoretical descriptions of this phenomenon fall into two 
groups. In the first, the Gross-Pitaevskii equations for the two 
condensate wavepackets are modified (either phenomenolog- 
ically 01, or on the basis of a method of approximating quan- 
tum field theory 01) to give an elastic scattering loss term. 
While these methods yield equations of motion which allow 
for depletion of the condensate wavef unctions, they do not 
include a description of the scattered atoms, and hence can- 
not describe the effects of bosonically stimulated loss. In the 
second class of treatments Vk[ the quantum field theory is 
linearized about the condensate, yielding equations of motion 
linear in the fluctuation operators. This method shows that the 
process is essentially one of four-wave mixing between the 
two condensate fields and pairs of quantized fluctuations — 
however, as a linearized theory, it can deal only with perturba- 
tively small amounts of scattering and cannot simultaneously 
account for depletion of the condensate. Both of these for- 
malisms are valid only in the limit of weak scattering, but fail 
for large scattered fractions, such as we treat in this paper. 

In this work we will show that a treatment in terms of a 
classical field with added quantum fluctuations is not only 
able to produce the scattering halo, but also predicts a hith- 
erto unobserved phenomenon, which we shall call quantum 
turbulence, in the resulting halo, which should be straightfor- 
wardly observable. This method is able to describe: a) the 



evolution (including large depletion) of the condensate wave- 
function as the scattering progresses; b) the full quantum cor- 
relations originating from the nonlinear amplification of the 
quantum vacuum fluctuations. 

The classical field method can be formulated as follows: 
i) The system is described by a c-number field amplitude 
^(x, t) which satisfies the Gross-Pitaevskii equation {I} and 
has a mode expansion 
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where aj (t) is the amplitude for the mode with wavevector kj 
and V = LxLyLz is the volume contained within the periodic 
boundaries in coordinate space. 

ii) The quantum fluctuations are introduced in the initial con- 
dition for ^(x, t), which is given by the sum ^ (x, t = 0) = 
{x.) -\-x (x) , where (x) and x (x) are respectively the real 
and virtual particle fields. We express the field of virtual par- 
ticles using X (x) = Zljli Xj (^kj • x) /VV, where the 
amplitude in each mode is Gaussian with the properties 
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The mean value of the total virtual population is thus M/2. 

The fundamental basis for the classical field method is the 
representation of the system of many bosons by means of a 
Wigner function description, which for l arge occupations is 

equivalent to this form lili[i3[iil[il[ii[i3[iE[il[il 

flii fl9[ [201 (the truncated Wigner approach). Heuristically, 
the introduction of virtual particles via Q can be viewed as 
adding half of one particle per mode, corresponding to the 
zero-point occupation of the ground state of the harmonic 
oscillator which represents each mode. The choice of ini- 
tial condition ^ is correct where V^(x) = 0, but gener- 
ates a slightly heated and nonequilibrium condensate where 
?/^(x) 7^ EUl . This will make very little difference to 
the results of the calculations, which require that the vacuum 
be represented accurately — furthermore, in practice the differ- 
ence between the state thus represented and a pure condensate 
is very much less than the experimental uncertainty. 



The pseudopotential approximation used in Eq.O results 
from the eUmination of modes with momenta larger than a cer- 
tain cutoff. In the work of [sll this cutoff effectively leaves 
two sets of modes, centered around the mean momenta of each 
of the individual colliding wavepackets, and mutual scattering 
appears explicitly. In contrast, we use a single spherical mo- 
mentum space cutoff which is large enough to include all rele- 
vant momenta for the scattering process, but small enough for 
the pseudopotential approximation to retain its validity. This 
is possible since the essential requirements for the elimina- 
tion of the higher modes are that these modes should have 
no occupation, and that the wavelength at which the cutoff is 
made should be larger than the range of the interatomic po- 
tential. The field is evolved using the Gross-Pitaevskii equa- 
tion within a projection formalism similar to that introduced 
inQ. 

For this problem, we write the Gross-Pitaevskii equation 
for each mode as 



-( 

2m 



mV 



^a*asat5jr,st, (4) 



10 

-10 

10 

-10 

I 10 

-10 



1.5 



1.0 



+ 



0.5 



where the function Sjr^st is unity for + — kg — = 
and zero otherwise. The projection is implemented by the 
requirement that all of |k^ | , |k^ | , Ik^ | , |kt | be less than i^max, 
the wavenumber which sets the momentum space cutoff. 

The nonlinearity present in the Gross-Pitaevskii equation 
gives rise to pairwise interactions of the type j-\-r ^ s-\-t, but 
onlyif the amplitudes a^, <^s and at are all non-zero. As noted 
in this mechanism can be viewed as a non-degenerate 
parametric oscillator, as found in nonlinear and quantum op- 
tics i22ll . In the classical field formulation, the noise terms 
Xj ensure that all modes do have non-zero occupation. Vac- 
uum modes, in which the occupation is initially virtual, can 
therefore develop a macroscopic population as a result of the 
interaction. Thus the scattering in our formalism arises natu- 
rally as a result of the classical field method. 

The usual derivation of the classical field formulation via 
the truncated Wigner function 

Hi El El 

requires that 

all relevant occupations be large, and does not immediately 
seem to be satisfied here since we are considering modes 
whose only occupation is the virtual population. However, 
for this situation, in which the most significant contributions 
come from terms in which two of the aj amplitudes in Q 
are macroscopically occupied, it is still possible to prove the 
negligibility of all terms with third-order derivatives in the 
Fokker-Planck equation which represents the exact equation 
of motion in the Wigner representation, and this is the condi- 
tion for the validity of the classical field method. This result 
will be published elsewhere. 

We consider collisions between two equally populated 
wavepackets derived from a single trapped condensate using a 
Bragg grating of very short duration compared to the overlap 
time of the packets 12^ 123]. Thus an appropriate initial state, 
t/j (x) of the real particles is the sum of two distinct momen- 
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FIG. 1: (a-f): Velocity mode populations on the planes Vz = 
(left) and Vx = (right) for the collision of packets generated from 
a six million atom condensate at t = (top), t — 0.5 ms (middle) 
and t — 2.0 ms (bottom). The spherical momentum cutoff is clearly 
visible in the upper plots due to the presence of quantum fluctuations, 
(g-h): Mode populations at t = 2.0 ms for an identical collision 
excluding vacuum noise. 



tum wavepackets with the same spatial envelope: 

V;(x) = -i^V^o (x) [e'- 



(5) 



where V^o (x) is the stationary solution of the Gross-Pitaevskii 
equation for the Nq atom condensate in a trap and k± = 
{±mvc/h, 0, 0) are the wavepackets' normalized amplitudes 
and wavevectors. For this paper we choose the case of a ^^Na 
condensate in an axially symmetric trap with o;^ = 27r x 80 
Hz and cj^ = 27r x 20 Hz, as used in fs*], but with a maxi- 
mum of 6 X 10^ atoms, compared to the maximum number 
3 X 10^ used in the experiment. Immediately following the 
application of the Bragg pulses (t = 0) the confining potential 
is removed, so that the evolution takes place in free space. 

The velocity space population distribution for the collision 
of packets from an initial condensate of 6 x 10^ atoms with 
Vc = 10 mm/s is shown in Fig. [ij a-f) for a sequence of times. 
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This figure shows: 

i) The original wavepackets broadening and changing shape 
from circular in projection to crescent shaped. 

ii) The generation of a circular feature — the scattering halo — 
centered at the system center of mass velocity with a radius in 
velocity space approximately equal to (but larger than) Vc. 

The scattering halo is characterized by patches of high 
mode population, the phase grains, whose size and aspect ra- 
tio are very similar to those of the original condensate packets 
(in velocity space). This characteristic size can be understood 
in terms of the parametric amplification of vacuum fluctua- 
tions described by Eq.(|4|). Pairs of vacuum modes j and r 
at each end of a diameter of the energy conservation sphere 
are selected preferentially for growth (or loss) if their phases 
are appropriate. The sum in Eq.(|4|) convolves the noise field 
near mode r with the initial condensate packets, and if mode 
r grows preferentially, it drives the growth of all modes near j 
within a volume of the size of the condensate velocity wave- 
function. 

In the weakly populated regions separating the phase 
grains, analysis of our data reveals a large number of vortices, 
which is indicative of a turbulent velocity field and arises di- 
rectly from the parametric amplification of the vacuum fluc- 
tuations. This justifies our identification of the behaviour as 
quantum turbulence. (An analogous amplification of electro- 
magnetic vacuum fluctuations was achieved in 1982 using a 
Josephson junction 12^1 .) 

In contrast, a simulation performed without the vacuum 
fluctuations generates a simple modulational instability, as 
shown in Fig.[Ilg-h). The fluctuations are equivalent to 
3 X 10^ additional particles, and thus their inclusion repre- 
sents a significantly altered initial condition from that given 
by omitting them altogether. 

As the halo generation process is analogous to parametric 
amplification, it is expected that modes of opposite velocity 
in the halo will display correlations, and that modes of simi- 
lar velocity will similarly display correlations due to the finite 
size of the phase grains. Calculating these correlation func- 
tions should properly be performed on an ensemble of sim- 
ulations. However, for this situation which, away from the 
condensate packets, displays a high degree of spherical ho- 
mogeneity in velocity space, we can use a single simulation 
and average over modes in velocity space which have similar 
speeds and which are not occupied by the condensate packets. 
To this end we define the data collection region, consisting 
of those modes whose polar angles, (^j, lie between 7r/4 and 
37r/4 to the positive Vx axis. Clearly such a region would also 
be needed experimentally. Within this region then, we define 
the correlation functions 
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M{v) = {aiaj)^^^^^, (8) 
as the average non-condensate mode population (excluding 
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FIG. 2: Correlation functions for the six million atom simulation 
averaged over the data collection region as described in the text: (a) 
average mode population; (b) pair creation correlation; (c,d) ampli- 
tude autocorrelation. 

virtual particles), the amplitude autocorrelation function and 
the pair creation correlation function respectively. The corre- 
lation functions for the six million atom simulation are shown 
through time in Fig.|2| 

From Fig.|3a), we observe that population builds up ini- 
tially near v = 10 mm/s, peaking ait = 0.52 ms for v = 9.3 
mm/s. The subsequent redistribution of population within the 
halo modes occurs via inter-halo scattering events and can be 
considered as a thermalisation process. 

Fig.|2jb) shows the pair creation correlation for the paired 
modes, = — Vj and shows that a definite, although small 
10%), correlation develops, peaking at t = 0.30 ms and 
then decaying as further scattering events occur. 

Fig.|3c,d) show the velocity space amplitude autocorrela- 
tion functions atv = 9.3 mm/s. Fitting a Gaussian in Svx,y,z 
to obtain the FWHM correlation lengths, we find that from the 
time that population is established (at about t = 0.05 ms) the 
ratio of the 6vx and Svy correlation lengths to the Svz length 
remain close to 4. The individual lengths initially grow and 
peak at t = 0.25 ms (where they are 1.2 mm/s in the Svx and 
Svy directions, and 0.3 mm/s in the Svz direction) and then 
decay slowly. These can be compared to the FWHM of the 
original velocity space condensate wavefunction, which is 0.7 
mm/s in the Vx and Vy directions and 0.18 mm/s in the Vz di- 
rection. This correlation behavior is a quantitative measure of 
the size of the phase grains, and we estimate the average num- 
ber of atoms in a coherent patch to be 150 at t = 0.25 ms for 
V = 9.3 mm/s. 

In Fig.|3lwe examine the time evolution of the real particle 
population in the data collection region, defined using 
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FIG. 3: Normalised data collection region populations for the classi- 
cal field simulations (sold lines) and complex scattering length sim- 
ulations (dashed lines). Initial condensate numbers are (lowest to 
highest) (1, 3, 6) x 10^, all with Vc = 10 mm/s. Subplot (b) shows 
the early times of (a). 

for varying initial condensate number. For larger condensate 
numbers a larger fraction is scattered, and the scattering is 
completed earlier. The decrease in population at later times 
is due to interhalo scattering events redistributing population 
outside the data collection region. 

For comparison, we plot the fraction of condensate popu- 
lation lost into the data collection region (assuming isotropic 
loss) during an identical collision using the complex scattering 
length method of |4]. The most significant difference using 
our method is a considerably larger out-scattered fraction, be- 
cause Bosonic stimulation is included. Additionally, the scat- 
tering model of 01 is essentially based on Fermi's "golden 
rule", which results in exact kinetic energy conservation at all 
times, and gives an initial linear growth of the out-scattered 
fraction. However, our method shows an initial quadratic 
growth rate, because the timescale over which one must av- 
erage to derive the kinetic energy conservation implicit in the 
golden rule (say to an accuracy of about 10%) in this case is 
about 0.5s. As a result of both of these effects, the rate of 
population loss, instead of being greatest at t = 0, where the 
wavepackets are maximally overlapped, is greatest at about 
t = 0.3 -0.5ms. 

The major loss mechanism in a BEC is that of three body 
recombination |27|,|28ji], but this can have very little effect on 
the results. Using the experimentally obtained loss rate, we 
estimate that in the time taken (^ 2ms) for an experiment, no 
more than 5 atoms would be lost from the system. Implement- 
ing a phenomenological three-body loss into our simulations 
shows no change in the dynamics until we increase the loss 
rate to more than 100 times the experimentally obtained value. 

Conclusions: This paper makes the first application of the 
classical field method to a realistic three dimensional prob- 
lem, and produces results which could not have been gener- 
ated by any other method in current use. Our detailed nu- 
merical results require only a single run of the simulation, 
avoidingtime-consuming ensemble averages. Previous calcu- 
lations 13, 0, H, 01 have been restricted to regimes in which 
Bosonic stimulation is not important, and have only been able 
to determine more limited information, whereas our method 
can handle very large (?^ 40%) outscattered fractions. 

The only significant restriction on our method is the com- 
puter memory required to represent the system as it expands 



in space. For the situation considered here this is not a prob- 
lem. Collisions with smaller out-scattered fractions would be 
more difficult to simulate because of the high ratio of virtual 
to real particle occupation, but the method would still be valid 
if sufficiently large ensembles of initial conditions were used, 
whereas for large out- scattered fractions, a single realization 
is sufficient for most measurable quantities. 

New qualitative features found are: i) The existence of 
quantum turbulence, which should be easily detectable ex- 
perimentally, using for example measurement of the density 
along a slice with the same techniques used to observe simi- 
lar slices through clouds containing vortex lattices 1131; 
ii) The suppression of the modulational instability which a 
simple Gross-Pitaevskii picture would predict; iii) The re- 
duction in correlations expected between scattered atoms of 
opposite momentum to a rather small value. 
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